Arithmetic means of speculative coverage are in Figure 1.

Figure 1: Arithmetic means of crop speculative coverage.

coverage.lm1 <- lm(speculative.pct.coverage   ~ as.factor(block) + species , data = coverage, subset = (treatment != "N"))
resid_panel(coverage.lm1)

joint_tests(coverage.lm1)

##  model term df1 df2 F.ratio p.value
##  block        3  30   8.800  0.0002
##  species     10  30  12.567  <.0001

plot(emmeans(coverage.lm1 , "species"), comparisons = TRUE) + 
#    scale_y_discrete(limits = rev(levels(species))) +
  xlab("Speculative coverage (%)" ) +
  ylab("Cover crop")  + 
#  coord_flip() +
  theme(text = element_text(size=20)) +
  theme_classic()

Figure 2: Estimate speculative percent coverage of different cover crop species

coverage.emm <- emmeans(coverage.lm1 , "species") %>% as.data.frame()

anova <- aov(speculative.pct.coverage   ~ as.factor(block) + species , data = coverage, subset = (treatment != "N"))

# Tukey's test
tukey <- TukeyHSD(anova )

# compact letter display
cld <- multcompLetters4(anova, tukey)
# https://www.mathiasecology.com/code 
cld.df <- data.frame(multcompLetters(TukeyHSD(anova)$species[,4])$Letters) 

cld.df$species <- rownames(cld.df)

coverage.emm %>% left_join(cld.df, by = "species")

##  species           emmean   SE df lower.CL upper.CL
##  Arugula             3.12 3.86 30    -4.75     11.0
##  Broadleaf mustard  23.75 3.86 30    15.87     31.6
##  Brown mustard      29.38 3.86 30    21.50     37.3
##  Collard            37.50 3.86 30    29.62     45.4
##  Kale               12.50 3.86 30     4.62     20.4
##  Mustard            35.00 3.86 30    27.12     42.9
##  Oriental mustard   31.25 3.86 30    23.37     39.1
##  Radish              5.00 3.86 30    -2.88     12.9
##  Rapeseed            3.12 3.86 30    -4.75     11.0
##  Red clover          3.12 3.86 30    -4.75     11.0
##  Turnip             17.50 3.86 30     9.62     25.4
##  multcompLetters.TukeyHSD.anova..species...4...Letters
##  e                                                    
##  abcd                                                 
##  abc                                                  
##  a                                                    
##  cde                                                  
##  ab                                                   
##  abc                                                  
##  de                                                   
##  e                                                    
##  e                                                    
##  bcde                                                 
## 
## Results are averaged over the levels of: block 
## Confidence level used: 0.95

coverage %>%
  filter(treatment != "N") %>%
  dplyr::group_by(treatment, species2) %>%
    mutate(sp.coverage = mean(speculative.pct.coverage),
    sd = sd(speculative.pct.coverage),
    n = n(),
    se = sd / sqrt(n)) %>%
  left_join(cld.df, by = "species") %>%
  ggplot(aes(x = species, y = sp.coverage)) +
  geom_bar(stat = "identity", position = "dodge", fill = "#B31B1B") + 
  geom_errorbar( aes(x = species, ymin = sp.coverage - sd, ymax = sp.coverage + sd), width=0.4, colour="black", alpha=0.9, size=1) + 
  geom_text(aes(y = 59, x = species, label = multcompLetters.TukeyHSD.anova..species...4...Letters)) +
  ylab("Speculative coverage (%)" ) +
  xlab(" ") + 
   theme(legend.title= element_blank(),
         axis.text.x = element_text(face = "italic")) + 
  scale_x_discrete(labels= c("Arugula \n (Eruca vesicaria)" ,
                            "Broadleaf mustard \n (Brassica juncea)",
                            "Brown mustard \n (Brassica juncea)",
                            "Collard \n (Brassica oleracea)",
                            "Kale \n (Brassica oleracea)",
                            "Mustard \n (Sinapis alba)",
                            "Oriental mustard \n (Brassica juncea)",
                            "Radish \n (Raphanus raphanistrum subsp. sativus)",
                            "Rapeseed \n (Brassica napus)",
                            "Red clover \n (Trifolium pratense)",
                            "Turnip \n (Brassica rapa)"),
                   guide = guide_axis(angle = 40)) # + scale_x_discrete(guide = guide_axis(angle = 40))